Dirac fermion quantization on graphene edges: 
Isospin-orbit coupling, zero modes and spontaneous valley polarization 
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The paper addresses boundary electronic properties of graphene with a complex edge structure of 
the armchair/zigzag/armchair type. It is shown that the finite zigzag region supports edge bound 
states with discrete equidistant spectrum obtained from the Green's function of the continuum 
Dirac equation. The energy levels exhibit the coupling between the valley degree of freedom and 
the orbital quantum number, analogous to a spin-orbit interaction. The characteristic feature of the 
spectrum is the presence of a zero mode, the bound state of vanishing energy. It resides only in one 
of the graphene valleys, breaking spontaneously Kramers' symmetry of the edge states. This implies 
the spontaneous valley polarization characterized by the valley isospin ±1/2. The polarization is 
manifested by a zero-magnetic field anomaly in the local tunneling density of states, and is directly 
related to the local electric Hall conductivity. 

PACS numbers: 73.20.At,73.22.Gk,73.63.Bd 
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I. INTRODUCTION 



Due to the close connection between their topological 
and physical properties, two-dimensional (2D) electron 
systems have traditionally been in the focus of funda- 
mental research. From the practical side, device func- 
tionalities in the 2D geometry are of great importance 
for applications and particularly suitable for lateral elec- 
tronic architecture. The interest in these general aspects 
of 2D electron systems has recently revived in the light 
of the experimental success in isolating individual lay- 
ers of graphite, preserving the honeycomb crystal struc- 
ture^. Such a system - graphene - exhibits elemen- 
tary excitations behaving at low energies and long dis- 
tances as massless Dirac fermions 3 ' 4 . Due to its massless 
quasiparticles graphene stands out among other 2D elec- 
tron systems, which is probably most prominently mani- 
fested by the unconventional quantum Hall physics (e.g. 
Refs. the phenomenon of Klein tunneling^ 

and fermion bound states on extended defects such as 
graphene boundaries^^^^ii^^iiii&i 9 .^, to name 
a few. In particular, understanding boundary effects 
in clean and disordere d 21 ' 22 graphene and the need for 
their characterization are among the outstanding current 
challenges in the field, arising from potentially promis- 
ing electronic applications of graphene ribbon o 23 ' 24 and 
quantum dots 2 ^. 

One of the reasons why the boundary effects in 
graphene should matter was pointed out quite a time 
ago by Fujita et al [Ref. Using tight-binding cal- 

culations they predicted a new branch of quasiparticle 
states localized on the so-called "zigzag" edge. It is 
one of the most common types of the honeycomb lattice 
termination formed by two parallel crystal faces of the 
triangular sublattices of the honeycomb structure [see, 
Fig. UKa)]. The properties of the zigzag edge states are 
better understood when compared to the edge states in 
conventional 2D quantum Hall systems^ 6 -. Unlike the lat- 
ter, the zigzag edge states exist without any external 



magnetic field and any excitation gap in the 2D bulk. 
They are nonchiral: there is a Kramers' pair of counter- 
propagating modes originating from two nonequivalent 
nodal points of graphene's Brillouin zone [see, Figs. IHb) 
and (c)] . The zigzag edge states have essentially the same 
origin as the bound states of massless fermions on domain 
walls^ 7 -. Here the role of the domain wall is assumed by 
the out-of-plane rotation of the "sublattice" spin which 
in the continuum limit corresponds to the zigzag edge2£. 
Experimental evidence for the bound states on graphene 
edges comes from both tunnelin g 12 ' 13 and angle-resolved 
photocmission spectroscopies 1 ^. 

The present study is motivated by the observation that 
in experiments one has to deal with finite-length zigzag 
edges that represent a section of the graphene boundary 
sided usually by two armchair edge s 12 ' 13 . As the arm- 
chair sides do not support edge state s 10 ' 16 , one should 
generally expect quantization of the propagating modes 
in the finite zigzag section. This type of quantization 
is distinct from the size-quantization in zigzag graphene 
ribbons studied earlie r 10 ' 11 ' 15 ' 16 , because it can occur on 
an isolated zigzag boundary, which is the typical situa- 
tion in scanning tunneling experiment o 12 ' 13 . The conse- 
quences of such a quantization have not been studied pre- 
viously. In the present work, they are addressed within 
the Dirac fermion confinement model derived from the 
lattice structure shown in Fig. HJa). 

In our approach the time-reversal symmetry and 
Kramers' degeneracy of the zigzag edge states comes as a 
result of an effective isospin-orbit coupling. The isospin 
is introduced as a convenient formal representation for 
the two nonequivalent nodal points of graphene's Bril- 
louin zone. The rotations generated by the isospin leave 
the 2D Dirac equation invariant. We show that this con- 
tinuous symmetry is broken by the zigzag confinement, 
and the edge state spectrum explicitly depends on the 
confinement parameters controlling the isospin-orbit cou- 
pling. The quantization of the edge states is achieved by 
imposing effective boundary conditions at the ends of 
the zigzag edge [see, Fig.QJb)]. They cause the interval- 
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FIG. 1: (Color online) (a) Example of a finite-length zigzag 
edge sided by two armchair boundaries. A and B mark the 
sites of the two triangular sublattices. (b) Geometry of the 
continuum model for the system in panel (a): The zigzag edge 
at y — supports a Kramers' pair of counter-propagating 
edge states from K+ and K- valleys. They transform into 
each other via intervalley scattering caused by the armchair 
sides at x — ±L/2. (c) Schematic view of the spectrum near 
the nodal points K+ and K- of graphene's Brillouin zone^. 



II. THE BOUNDARY PROBLEM 

A. 2D massless fermions, chiral symmetry and 
isospin 

The two distinct nodal points (valleys) of graphene's 
Brillouin zone result in a pair of massless Weyl fermions 
whose wave functions, tp + and ^>_, satisfy the matrix 
equation: 



,H 



crp 







U(crp)U- 



•(1) 



It is assumed that the Hamiltonian H is diagonal in val- 
ley space (+, — ). The intra- valley Hamiltonians are ex- 
pressed in terms of the Pauli matrices cr XjVtZ acting on 
the functions 



(2) 



that have two components due to the bipartite lattice 
structure of graphene, with two sublattices denoted as A 
and B in Fig. [IJa); v and e are the Fermi velocity and 
energy with respect to the Fermi level, and the quasipar- 
ticle momentum p is confined to the plane of the system. 

We further assume that the intra- valley Hamiltonians 
are related to each other by the chiral symmetry: 



[/(crp)?/- 1 = -crp, 



(3) 



ley scattering connecting the incident and outgoing edge 
states, which models the armchair confinement. It turns 
out that the quantized spectrum contains a zero mode, 
i.e. the state with vanishing momentum and energy. Re- 
markably, it couples only to one of the isospin projec- 
tions, that is it exists only in one of the valleys, breaking 
spontaneously the Kramers' symmetry of the edge states. 
This leads to the spontaneous isospin (valley) polariza- 
tion with the total edge-state isospin ±1/2. This mech- 
anism of the valley polarization differs from the previous 
proposals^. We demonstrate that the spontaneous sym- 
metry breaking can be detected through the magnetic- 
field dependence of the tunneling density of states, and 
also find a direct relation between the isospin polarization 
and the local electric Hall conductivity. 

The subsequent sections give a complete account of our 
approach: In Sec. |TT] we formulate the boundary prob- 
lem for a finite zigzag edge and analyze it in terms of 
the discrete and continuous symmetries of the 2D Dirac 
fermions. The Green's function solution of the boundary 
problem and the spectrum of the quantized Dirac fermion 
edge states are discussed in Sec. Mil Section llVl addresses 
the valley polarization effects, both spontaneous and in- 
duced. The latter is the analogue of the quantum spin 
Hall polarization. Finally, section [V] describes the sig- 
natures of the valley polarization in observables, such as 
the tunneling density of states and the local electric Hall 
conductivity, and contains concluding discussion. 



where U is a unitary matrix. In this way we explicitly 
account for the generic property of nodal lattice quasi- 
particles known as fermion doubling: they come in pairs 
of opposite-chirality (Weyl) species that together obey 
the Dirac equation^. We note that in the 2D case the 
unitary transformation, Eq. ([3]) is always achieved by one 
of the a matrices. If, for instance, the system is located 
in the x, y plane [Fig. [T], we have 

crp = a x p x + GyPy, U = a z . (4) 

The discrete chiral symmetry, Eq. Q can be promoted 
to a continuous one. Let us introduce another set of 
the Pauli matrices ti,2,3, acting in the valley space, and 
consider the vector operator, 

1=2 ® O*, T2 <g> a z ,T3 ® (JQ^j , [Ik,Il] = iSklmlm, (5) 

whose components Ik (k — 1,2,3) formally satisfy the 
commutation relations of an angular momentum. It is 
easy to see that the Hamiltonian H is invariant under 
rotations generated by Ik'- 

e ie k i kHe -iO k i k = H H = VT 3 ®(crp), (6) 

where 8k is the rotation angle. This means that the 
original choice of the upper and lower components of 
* [Eq. flU] as being the "+" and "-" valley functions, 
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respectively, is not physically distinguished. One can 
rather treat them as the " up" and " down" states of the 
effective spin (isospin) 1/2. We will nevertheless keep 
the original notations rp± for the upper and lower com- 
ponents of 'J, interpreting them as the projections 

lk=(!±/ 3 )*, (7) 

where Iq = tq ® er is the 4x4 unit matrix (the direct 
product of the 2x2 unit matrices To and o~q). 

B. Boundary condition for the zigzag edge and 
broken isospin rotation symmetry 

The zigzag edge is a type of the honeycomb lattice 
termination where the outermost lattice sites all belong 
to one of the sublattices [Fig. [T{a)]. It does not couple 
the states from the K + and K- valleys^, due to which 
the continuum boundary condition for the zigzag edge 
can be obtained by rather simple reasoning^ To be con- 
crete let us assume that the outermost sites are all of 
the A type and the next (missing) atomic row would be 
of the B type, as in Fig. [Ha). O n the missing B row 
one can impose the hard-wall condition ipB±(x,0) = 0, 
while keeping if>A±(x, 0) arbitrary. In spinor notations 
[Eq. ©], this reads 

^±M) = o4l^±M), 1x = (0,0,1). (8) 

This boundary condition admits the generalization be- 
yond the hard-wall approximation. It is achieved by ro- 
tating the unit vector lj_ about the normal to the 
boundary (in Fig. QJb), iib|| y), which is consistent with 
the requirement for the normal component of the cur- 
rent to vanish at the edg o 18 i 30 . Moreover, the rotation 
can be made valley-dependent: lj_ — > 1±. Using the 4 
component spinors, we can therefore write: 

*(x,0) =M*(x,0), (9) 

, , T +T 3 T0—T3 , . 

AI = — - — (K)<rl + H — g>crl_, (10) 

4 = 1, (l±n B ) = 0. (11) 

Further restrictions on l± are imposed by the discrete 
symmetries of the problem. As the lattice prototype of 
our system has two identical sides [Fig. QJa)], our con- 
tinuum model should inherit spatial parity with respect 
to coordinate reflection along the edge, i.e. x — > —x. 
It is the symmetry of the Dirac equation |T]) since the 
coordinate reflection can be compensated by the spinor 
transformation: 

y p (x,y) = AV(-x,y), A = n®cr x , (12) 

simultaneously swapping both the valley and sublattice 
spinor components. However, the boundary condition, 



Eq. ([9]) does not share this symmetry because M and A 
do not commute 

AM A' 1 = T °~ T3 ® (a x l s + - <rj„+) + 

+ T ° 2 T3 ® ( a * l x- ~ (13) 

unless there is a relation between 1 + and 1_ such that 

l x+ = l x - = l x , l z+ = -l z - = l z , 1= (Z X ,0,Z 2 ).(14) 

These restrictions also make the zigzag boundary invari- 
ant under time-reversal operation ty T {x, y) = A^*(a;, y). 

We are now prepared to prove that the zigzag bound- 
ary condition, Eq. ([9]) violates the isospin rotation sym- 
metry. More specifically, we are talking about the non- 
trivial rotations generated by the I\ and I2 components 
of the isospin, Eq. ([5]) . Indeed, the matrix M (jTUJ) does 
not commute with lip: 

h,2 MI l,2 = T ° 2 T3 (~ a x lx + + < J zlz+) + 

+ ^^Vi-aJ^+aJzJ), (15) 
unless 1+ and 1 satisfy the conditions: 

lx+ — —l X -, lz+ = h.-- (16) 

These are incompatible with the requirements for the 
x-parity and time-reversal symmetry, Eq. (|14p . In sec- 
tion UII B\ we demonstrate that the broken isospin rotation 
symmetry implies an analogue of the spin-orbit coupling 
controlled by the components of the vector 1 in Eq. 

C. Parity-symmetric armchair edges 

We now turn to the boundary conditions at the arm- 
chair sides x = ±L/2. They should account for the val- 
ley and sublattice mixing specific to the armchair lat- 
tice termination^ and, at the same time, possess both 
the x-parity and time-reversal symmetry. The suitable 
boundary conditions can be written as^i 

*(±f,y)=A*(±f,y), (17) 

with the same off-diagonal matrix A as in Eq. (I12|) . They 
meet the requirement of the vanishing of the normal com- 
ponent of the Dirac current: 

j x (±±,y) = & (±±y)T 3 ®a x *(±±,y) (18) 
= *t (±f ,y)n ®<t x (t 3 $ a x )n ®a x ^ (±f ,y) 

= -&(±% ) y)-n®v x *{±%,v) = -jz (±|,y) =0, 

where we have switched to the creation ^ (x, y) and an- 
nihilation ty(x,y) operators. 

Importantly, the x-parity of the problem allows us to 
reduce the boundary conditions, Eq. (TIT)) to the usual 
symmetric boundary conditions: 

= *(-§>»)■ (!9) 
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To prove this we first notice that the original function 
fy(x,y) and the transformed one ^ p (x,y) [Eq. (fT2"j) ] cor- 
respond to the same solution of Eqs. H}, @ an d (fl7|) , 
and, therefore, must coincide: ^(x,y) = A^>(~x,y). In 
particular, at x — ±L/2 we have 



*(±f,y)=A* ( T %,y) 



(20) 



Comparison with Eq. (fT7|) yields Eq. (fTO")) . For the ac- 
tual calculations, we will use the symmetric boundary 
conditions modulated by a magnetic phase <j>: 



v 2 ' 



y) = * 
E m = 



L 
2 • 



y) exp(2 7 r l 0), (21) 
(h/eL)4>. (22) 



III. 



In this way we account for a weak magnetic field per- 
pendicular to the plane x,y. If its vector potential is 
chosen to be parallel to the zigzag edge, A(y)||x and to 
vanish at y — > oo, then at y — the phase <p exactly 
equals to the flux through the strip in units of ch/e. For 
weak magnetic fields, the spatial variation of cf> with the 
coordinate y can be neglected, while its adiabatic varia- 
tion with time implies an electric field along x given by 
Eq. j22]). 



DIRAC FERMION EDGE STATES ON 
FINITE ZIGZAG EDGES 

A. Green's function of the system 



To study the spectral properties of zigzag graphcne 
edges it is convenient to use the Green's function ap- 
proach. The specifics of its implementation to boundary 
problems i n grap hene is still scarcely covered in literature 
(e.g. Refs. fl7il33 ). Below we describe in some detail the 
main calculation steps leading to the final result given by 
Eqs. flU - AMD- 

We begin by introducing the retarded Green's func- 
tion G R (rt,r't') as a 4 x 4 matrix in space of the valley 
(isospin) and the sublattice degrees of freedom whose ma- 
trix elements are given by 



G R (rt,r't') 



0{t - 1') 
ih 



x U a {vt)^\{v't') + ^(rY)V> Q (ri)) , 



(23) 



where the brackets (...) denote averaging with the equili- 
birum statistical operator and the indices a and (3 inde- 
pendently run through all possible combinations of the 
isospin and sublattice indices: a, (3 = A+, A_, _B + , £?_. 
As the zigzag edge [Eq. (|9])] possesses the isospin rota- 
tion symmetry generated by I3 (i.e. does not couple the 
valleys), G R can be decomposed into the direct product: 

G R (rt, rV) = \ J2 fa> + tt 3 ) ® G R (rt, v't') (24) 

1 r=±l 



where r = ±1 labels the valleys (i.e. the two isospin 
projections) and 



G R (rt,r't') = 



G AAlT (vt,v't') G AB[T (rt,r't') 
G BA \ T (rt,r>t') G BB]T {vt,v't') 



(25) 



is the matrix Green's function in sublattice space. Its 
time Fourier transform satisfies the equation 



(e<r - vT*p)G R (r, r') = <j 6(v - r'). 



(26) 



In terms of G R (r,r r ) the boundary conditions, Eqs. © 
and (|21[) of the previous section, read 



G?=(<rlr)G*\ y=0 , l T = l±, (27) 
G R \ x=L/2 = G R \ x= _ L/2 exp(27r#). (28) 

The solution to Eq. (|2"6")) can be sought in the form 



Gf(r,r') = U + — <xp x 



E 

TlgZ 



G AA \ T k n (y,y') 

G BB \ Tkn (y,y') 



(29) 



where the diagonal matrix elements are the Green's func- 
tions on sublattices A, B . They are expanded in plane 
waves e iknX with the wave number 

k n = {2ir/L)(n + (f>), n G Z (0, ±1, ...), (30) 



given by the boundary condition, Eq. (|28[) . For 
G AA:BB \ T k n {y, y') one has the ordinary differential equa- 
tion, 



(dy - ql)G A A,BB\Tk„(y,y') 



h 2 v 2 



S(y-y'), (31) 



and the boundary conditions following from Eq. (1271) 
"re(l - l ZT ) 



dyG AA \ T k„ 



9yG BB \ T k n 



HvL 



-Te(l+l ZT ) 



HvL 



G 



AA\rk„ 



y=o 



G 



BB\rk r , 



(32) 



, (33) 



y=o 



where q n = y/ k 2 — e 2 /h 2 v 2 . We seek the solution (finite 
at y — > 00) in the form 



G 



AA,BB\t 



kn (y,y')=C A , B (y'>- qny 



2h 2 v 2 q n 



-ln\y-y'\ 



where the first term is the solution of the homogeneous 
equation (f3Tj) and the second one is the Green's function 
of the unbounded system. The coefficients C A>B are ob- 
tained from Eqs. ([32]) and (|33|) with the following results: 



e -in(y+y') _ Q-in\y-y'\ 



G AAirkAy,y') = w ^ n 

| (1 + lzr)(q n + kg) - Tel XT /Hv ^_ qn[y+y r ) 
2(e — hvTk n l XT + iO) 



(34) 
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GBB\rk„(y,v') 

(1 

+ 



f e -Qn(,V+y') _ e —Qn\v-v' 



2fr 2 v 2 q n 

(1 ~ l Z r)(qn - k n ) + Tel XT /hv qn ( y+y ') 

2(e — twTk n l XT + iO) 



(35) 



The results of the above calculations can be summa- 
rized in the expression for the full matrix Green's func- 
tion: 



T + TT 3 



VT 

cq H (rp 



G*(r,r')= X! 

r=±l,ne2 

x (G* fc „ (y, y')* + G a rkn (y, y')a z ) ^ - ■ (36) 

We introduce the symmetric G s Tkn (y,y') — {GAA\rk n + 
G , ss|rfc„)/2 and asymmetric G* kn {y,y') = (G AA \ rkn - 
G BB\rk n ) /2 sublattice functions given explicitly by 



2n z v z q n V 

j On ~t~ k n l ZT 

2(e — hvrk n l XT + iO) 



-Qn(V+v') 



G a (y y') ^ n ^ n ^ ZT T~el XT / hv ^ 

r " ' 2(e — hvTk n l xr + iO) 



-in(v+v') 



(37) 



(38) 



In the equations above the denominators vanish at e = 
h,VTk n l XT . To identify this as a pole, we should make sure 
that the nominators remain finite as e — > hvrk n l XT . In 
this limit, the Green's function (l36l) behaves as 



G fl (r,r') « -— ^ {t + tt 3 ) ® (<t q + ct1 t ) 

r=±l,n£Z 

x Sll_ 9v e-l*«'-l(^')^(-»'), (39) 

e — nvTk n l XT + iO 

and we can see that the pole exists only if the unit step 
function Q(k n l ZT ) is not zero: 



hvTk n l xr , 



k n lzr > 0. 



(40) 



This is the spectrum of the states, decaying exponentially 
from the edge y = and propagating along x. 



B. Edge-state spectrum, isospin-orbit coupling and 
zero modes 



Let us analyze the edge-state spectrum, Eq. (|40"]) in 
some more detail. With the requirements of the x-parity 
and time-reversal symmetry [see, Eq. (|14|) ] and for k n 
given by Eq. ([ST)]) , we have 



e T . n = sgn(l x )AT(n + 
A = hv\l x \/L, 



T(n + <f>)l z >0, (41) 
n = 0,±l... (42) 



It is equidistant with the level spacing A and particle- 
hole asymmetric because of the restriction (n + (j))rl z > 



°+,n 



(^>0 



£ -,n 



1 2 3 4, n 

NT " 



-4 -3 -2 -1 




-n 



(()< ^-,n1 

1 2 3 4 vfl 



-4-3-2-1 |0 




n 



FIG. 2: Edge states in + and — valleys, Eq. flTj) for (a) pos- 
itive and (b) negative magnetic flux <f). We assume the zigzag 
confinement parameters, l x < and Z z > 0, so that the edge 
states exist below the Fermi level as inferred from tunneling 



spectroscopic measurements 



The external magnetic flux 



shifts the levels in + and — valleys in the opposite directions 
such that the zero mode n = (filled circle) occurs only in 
one of the valleys: + one for <j> > and — one for cj> < 0. 
The valley-dependent zero mode violates Kramers' symmetry 
of the edge-state spectrum for arbitrary small cf>. 



[see, Fig. |2J. The phase <j> results in the homogeneous 
shift of the levels. Let us consider \cj>\ <C 1 and neglect 
the shift in all of the states except the zero mode n = 0: 



£r,n = sgn(l x )Arn, rnl z > 0, n = ±1, ... (43) 
e Ti0 = sgn(/ a; )Ar0, T<f>l z > 0, n = 0. (44) 



We see that the states with n = ±1, ... exhibit Kramers' 
symmetry under r, n — > — r, — n resulting from the cou- 
pling between the valley (isospin) degree of freedom t 
and the orbital quantum number n. The isospin-orbit 
coupling originates from the broken isospin rotation sym- 
metry discussed in Sec. \II B\ The coupling constants are 
given by the parameters l x and l z of the zigzag confine- 
ment. For the hard- wall zigzag edge (l x = 0), we find 
the degenerate zero-energy state e r>n = 0. This is in 
agreement with the tight-binding calculations for zigzag 
graphene ribbons (e.g. Refs. llOillT ) if their results are 
extrapolated to the case of the infinite width when the 
edges become isolated. 



The zero mode, Eq. (j44[) stands out because it is due 
to the coupling between the isospin and the electromag- 
netically induced momentum fco = (2w/L)<f>. This mode 
breaks the Kramers' symmetry of the edge-state spectrum 
since it exists only for one of the isospin projections 
t — sgn(a^ z ), i.e. only in one of the valleys. In other 
words, there is a valley polarization effect. It is studied 
quantitatively in the next section. 
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IV. VALLEY POLARIZATION 
A. Spontaneous polarization 

To quantitatively characterize the valley polarization 
effect we introduce the local isospin polarization: 



P(e,r) 



ilmTr/ 3 G fl (r,r) 

7T 



= -^Z J2 ~I™G s Tkn (y,y), (45) 



T=±l,n& 



where Im denotes the imaginary part, the trace Tr of the 
Green's function ([55]) is taken in r <g> a space, and the 
function G s Tk ^(y,y) is given by Eq. (j37|) of the previous 
section. As we are interested in the edge isospin polar- 
ization, we relevant contribution to Im.G s Tk (y,y) comes 
from the pole in Eq. (|37)) : 



(46) 



x S(e 



0- 



Note that for the zero mode the step function Q(c/)tI z ) 
indicates the breaking of the Kramers' symmetry. 

Next we calculate the zero-temperature ground-state 
isospin density localized at the edge as 



i e {y) = / de Pe (e iy ) = (47) 

J — OG 

= -j d v E ^~ 2lknl ' lv Q(Kri z )e(-k n Ti x )(48) 



d v E ^^ k ^ye(k n rl z ). (49) 



To obtain the last formula we used the identity 
Q(x)Q(y) = Q(x)Q(xy). The summations over r = ±1 
and n can be done exactly: 



ie(y) 



N 
2L 



N n 



sgn <p e 
N = G(-l x l z )sgnl z , 



d y J2 sgn(n + 0)e- |n+ ^ |y/A 

n— — oo 

<t>\ y /x _ 2sinh(0y/A) 



A = 



ef/ A 
L 



Air . 



(50) 

, (51) 
(52) 



In Eq. ([51]) the first term, nonanalytic in <f>, is due to the 
zero edge mode n — 0. Its penetration length depends 
on the flux (f> and diverges at <j> — > 0. The second term 
accounts for the rest of the edge states n = ±1, ... It is 
an analytic function of <fr. The edge-state penetration 
length is measured in units of A given in Eq. (|5l2j) . and 
the flux is confined to a one-period interval chosen as 
-1/2 < <j) < 1/2. 

We note that depending on the boundary parameters 
the factor N (jB"2"]) takes integer values or ±1. The case 



N = corresponds to edge states above the Fermi level 
e = 0, for which the zero-temperature occupation number 
0(— l x lz) — 0. In what follows, we focus on the opposite 
situation, i.e. the edge states below the Fermi level and 



N = sgnl z = ±l, l x l z <0, 



(53) 



which is supported by the tunneling spectroscop y 12 ! 13 . 

Finally, we obtain the total isospin carried by the edge 
states as 



h. = L 



dyi e {y) = (^Y~ ~ <f>j ^gnl z 



(54) 



The zero mode results in the discontinuity at <f> = 0, due 
to which in the limit <fi — > the total isospin remains 
finite (half-integer): 



•sgn(0Z z ), 0^0. 



(55) 



This implies that the ground state does not share the 
time-reversal symmetry of the original equations \26\ ) 
H28\) in the limit <f> — ► 0. In this sense, the zero mode vi- 
olates the time-reversal and Kramers' symmetries spon- 
taneously, with the resulting spontaneous valley polariza- 
tion. 



B. Valley Hall polarization 

The accumulated isospin, Eq. (|54[) contains a linear 
term oc <j). It comes from the Kramers' degenerate edge 
states with n = ±1, ... [see, Eq. (|43|) ] . Such a prop- 
erty of Kramers' degenerate edge states was first no- 
ticed in the theory of quantum spin Hall systems (e.g. 
Refs. [33.34.35,36). The recent interest in these sys- 
tems is motivated by the principal possibility to real- 
ize a time-reversal invariant integer quantum Hall state 
in which the spin Hall conductance is quantized. From 
Eq. ([54]) it is possible to derive the analogue of the quan- 
tum spin Hall conductance. Let us calculate the isospin 
current as the rate of adiabatic change of the isospin: 
I e = 4>dl e /d(j) = —(eE x L/h)dI e /d<fi, which assumes 
4> 7^ and Eq. (|22|) . The derivative I e gives the trans- 
verse isospin flow in response to the voltage drop E X L 
along the edge: 



Ie — GmE x L 



Gm = -rsgnl z 
h 



(56) 



where the quantum isospin Hall conductance Gm takes 
the universal values ±e//i. As the zigzag-terminated 
graphene supports the edge states without any excita- 
tion gap in the bulk, the conductance (|56|) is hardly the 
signature of any bulk topological order—. We rather in- 
terpret it as the measure of the valley polarization rate 
at the edges. 
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FIG. 3: Tunneling density of states, Eq. (|60|l vs. magnetic 
flux (in units of ch/e) at various temperatures; y — 0.5A, 
i/o = 1/ALX. 



FIG. 4: (Color online) Tunneling density of states, Eq. (|60|l vs. 
temperature (in units of level spacing A) for small (<f> = 0.01) 
and large (4> = 0.5) flux values; y — 0.5A, vq — 1/ALX. 



V. SIGNATURES OF THE VALLEY 
POLARIZATION IN OBSERVABLES 



A. Tunneling density of states 



The presence of the valley polarization can be inferred 
from the magnetic field dependence of the zero mode. 
One of the possibilities is to measure the local tunneling 
conductance in the presence of a weak magnetic. At zero 
bias and the hnite temperature T, the tunneling conduc- 
tance is proportional to the tunneling density of states 



T7{T,</>) 



de 



df(e,T) 
de 



v(e,4>), (57) 



which is the convolution of the local spectral density of 
states, i/(e, 4>) and the energy derivative of the Fermi dis- 
tribution function, /(e,T). The local spectral density of 
states is obtained from the Green' function, Eq. ([36]) as 



f(e,r) 



-ImTrG fl (r,r) = 
~ ]T lmGl kn {y,y). (58) 



r=±l,7i63 



The edge-state contribution to lm.G s Tk (y, y) comes from 
the pole in Eq. ([371): 



v e {e,y) = -\d y e- 2 ^yQ(k n rl z )S^-^,n)-(59) 



From Eqs. ((57|) and (|59|) one can obtain the edge-state 
contribution to the tunneling density of states as 



5v(T,<t>) 



\<f>\v/>~ 



oo r 

E 



4iTA\ cosh 2^|A 

(n + ^e^y^ (n- 0)e^A 



(60) 



cosh' 



cosh 2 



where A is the level spacing given by Eq. ([42]) . 

Figure[3]shows that the flux dependence of Sv is nonan- 
alytic, indicating the spontaneous valley polarization at 
\4>\ — > 0. The nonanalyticity is present in a wide range of 
temperatures. The reason is that the zero-mode term al- 
ways dominates the flux dependence near = because 
it is linear in while the rest of the sum varies as 4> 2 . As 
demonstrated in Fig. [H for small (j) = 0.01 (red curve) 
the zero-mode also dominates the low-temperature be- 
havior of 8i>, showing a 1/T increase when T becomes 
much smaller than the level spacing A. This feature is 
due to the fact that for iji < 1 the energy of the zero 
mode oc </>A A. In contrast, for the rest of the sum 
in Eq. (f60|) the relevant energy scale is set by the level 
spacing A. 



B. Local electric Hall conductivity 

Although the zero-mode behavior in the tunneling den- 
sity of states signals the valley polarization effect, this 
observable does not provide the direct access to the ac- 
cumulated isospin. Here we intend to show that the ac- 
cumulated isospin is directly related to the local electric 
Hall conductivity. 

As the first step, we use Eq. (JSHJ) to calculate the zero- 
temperature ground-state charge density localized at the 
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FIG. 5: Hall conductivity in units of e 2 /h, Eq. {67J vs. (a) 
position and (b) dimensionless magnetic flux. 
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e <& e A A A Zigzag edge 





edge: 



Pe(y) 



deu e (e,y) 



(61) 



-d v Y, e- 2 ^y&(k n rl z )e(-k n Tl x ) (62) 
e6(-U z ) 



r,n£Z 



-2|/c„Z z |t/o 



Q(k n rl z ). (63) 



Again the summations over r = ±1 and n can be done 
explicitly: 



e9(-U z ) 

p e (y) = - T d, 



OO 

„ Y ^ n+ ^ v/X = (64) 

71— — OO 



eO{-lJ z ) 
L 



e -\<t>\y/X + 



2 cosh(</> y/X) 
ef/ A - 1 



,(65) 



with |0| < 1/2. The 4> dependence of Eq. (frJS*)) allows 
us to take the adiabatic time derivative, p e = ifid^pe = 
— {eE x L/h)d,f,p e , and obtain the following continuity 
equation: 



Pe *^y3y> Jy — ^yxE x , 



(66) 



where j y is the Hall current density induced by the 
transverse electric field E x , and $ yx is the local position- 
dependent Hall conductivity given by 



<h)x{y,4>) 



4e 2 L 



y / dy' i e (y',(f>). 



(67) 



It is expressed in terms of the edge isospin density given 
by Eq. (f5Tj) . The existence of the electric current density, 
j y normal to the system's boundary is consistent with 
the charge conservation because at the edge y = the 
conductivity <; yx {0, is zero [see, also, Fig.JSJa)]. It also 
vanishes far away from the edge: i yx (y — > oc, 4>) — > 0, so 
that the total edge charge is conserved: J °° dy p e = 0. 

At distances smaller than the characteristic penetra- 
tion length, y <C A, the Hall conductivity is simply pro- 
portional to the total isospin carried by the edge states: 

*.(»,*) « — —jUt) = -yj I — -<t>) ,(68) 



FIG. 6: (Color online) Other realizations of a finite-length 
zigzag boundary between two armchair edges. 



showing the same nonanalytic flux dependence as I e (4>) 
in Eq. §4$ [see, Fig.^b)]. 

In conclusion, we discuss the applicability of the re- 
sults of this paper. First of all, the lattice prototype 
of our continuous model [Fig. [lja)] is only one of many 
possible realizations of a finite-length zigzag boundary. 
It is nevertheless clear that for the honeycomb lattice an 
armchair/zigzag/armchair edge structure is rather typi- 
cal [see, e.g. Fig. [6]. Independently of its concrete real- 
ization, the edge states must experience multiple inter- 
valley backscattering from the two armchair regions, re- 
sulting in the bound states. The experimental estimate^ 
of the typical length of zigzag edges is of order of 10 ran. 
This is large enough for the applicability of our contin- 
uum model and, on the other hand, is shorter than the 
typical mean free path in graphene, which is required 
for the ballistic quantization. For samples with longer 
edges, multiple electron scattering due to boundary and 
bulk disorder may come into play, as revealed by recent 
numerical studie a 21 i 22 . The quantization effects studied 
in this paper are characteristic to isolated zigzag edges 
as opposed to the size-quantization in zigzag graphene 
ribbons. Because the control over graphene edges is still 
a serious experimental issue, it seems easier to obtain 
graphene samples with isolated zigzag edges rather than 
to produce zigzag-terminated ribbons. In this sense, lo- 
cal tunneling spectroscopy is currently the most adequate 
tool for investigating the edge states in graphene. As for 
the results on the local electric Hall conductivity (|67[) . 
their verification may present a challenging experimen- 
tal task. It should however be achievable with increasing 
control over the boundary effects in graphene. 
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